load timeseries;
%N=2^nextpow2(length(timehist));
N=length(timehist);

timehist=timehist/(1e6*60*60*24*365);
timehistint=linspace(0,timehist(end),N); dt=timehistint(2)-timehistint(1);

surfacehist=interp1(timehist,surfacehist,timehistint); surfacehist(1)=surfacehist(2);

for i=1:length(depthrange)
    qhist(:,i)=interp1(timehist,qhist(:,i),timehistint); qhist(1,i)=qhist(2,i);
    temphist(:,i)=interp1(timehist,temphist(:,i),timehistint); temphist(1,i)=temphist(2,i);
    avhist=interp1(timehist,avhist,timehistint); avhist(1)=avhist(2);
    %plot(timehist,qhist(:,i)); pause(1);
end

rate=diff(surfacehist);

thermaldepth=5; lowcut=1/100;

figure(1); clf; ax=plotyy(timehistint,surfacehist,timehistint,temphist(:,thermaldepth)'); set(ax(1),'ydir','rev'); linkaxes(ax,'x');
figure(2); clf; ax=plotyy(timehistint,bandpass(surfacehist,lowcut,inf,dt),timehistint,temphist(:,thermaldepth)'); set(ax(1),'ydir','rev'); linkaxes(ax,'x');
figure(3); clf; ax=plotyy(timehistint,bandpass(surfacehist,lowcut,inf,dt),timehistint,bandpass(temphist(:,thermaldepth)',lowcut,inf,dt)); set(ax(1),'ydir','rev'); linkaxes(ax,'x');
figure(4); clf; hold on; plot(temphist(:,thermaldepth)'-mean(temphist(:,thermaldepth)),surfacehist-mean(surfacehist),'b.'); plot(temphist(:,thermaldepth)'-mean(temphist(:,thermaldepth)),bandpass(surfacehist,lowcut,inf,dt),'r.'); plot(bandpass(temphist(:,thermaldepth)',lowcut,inf,dt),bandpass(surfacehist,lowcut,inf,dt),'g.');